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We develop a numerical method for approximating the surface modes of sphere- like nanoparticles 
in the quasi-static limit, based on an expansion of (the angular part of) the potentials into spherical 
harmonics. Comparisons of the results obtained in this manner with exact solutions and with a 
perturbation ansatz prove that the scheme is accurate if the shape deviations from a sphere are 
not too large. The method allows fast calculations for large numbers of particles, and thus to 
obtain statistics for nanoparticles with random shape fluctuations. As an application we present 
some statistics for the distribution of resonances, polariziabilities, and dipole axes for particles with 
random perturbations. 



PACS numbers: 41.20.Cv. 73.20.Mf 



I. INTRODUCTION 

The excitation of surface plasmons can cause strong in- 
teraction between light and metallic nanoparticles. These 
plasmons are hybrid modes of the electromagnetic field 
and the electron gas and are confined to the surface of 
the particle. They give rise to an enhancement of the 
incident field by several orders of magnitude [1-3] . This 
enhancement enables a variety of applications ranging 
from the well-established surface-enhanced Raman spec- 
troscopy (SERS), which allows the detection of even a 
single molecule [4, 5], to the emerging field of plasmon- 
ics [6, 7], which has led to prototypes of plasmonic wave- 
guides which effectuate optical energy transfer below the 
diffraction limit [6, 8, 9], 

A simple realization of a plasmonic waveguide is a 
chain of metallic spheres. A surface plasmon mode, typi- 
cally a dipole mode, of the first sphere of the chain is ex- 
cited and the scattered field of this first particle excites 
a surface mode in a sphere nearby and so the excita- 
tion can travel through the chain. There are two crucial 
points for an efficient transport: The spatial structure 
of the scattered field in the region of the neighboring 
sphere must allow for an efficient excitation of the fa- 
vored mode, and the overlap of the resonances of the 
bordering spheres has to be big enough. Since any real- 
ization of a sphere will deviate from an ideal one, thus 
introducing random fluctuations, it is important to esti- 
mate the typical magnitude of such deviations which still 
allow for an efficient transport. Therefore a simple and 
efficient numerical method for approximating the surface 
modes of the sphere-likes particles is needed. 

There are many numerical methods for the determina- 
tion of the electromagnetic field in the present of nano- 
sized scatterers, like the finite difference time domain 
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approach (FDTD), or so called semi-analytical methods 
based on expansions into special function systems like 
the multiple multipole method (MMP), or the discrete 
dipole approximation (DDA), sec, e.g., [10, 11] for re- 
views. Essentially, all methods able to calculate the fields 
also allow to determine the surface modes. For example, 
in [12] the DDA is used for the determination of the sur- 
face modes of nanoparticles, and in [13, 14] a boundary 
integral approach is proposed, which focuses on the sur- 
face modes, and has been used in [15] to determine the 
surface modes of single and coupled spheres, cylinders 
and cube-like nanoparticles. 

Here we use a semi-analytical approach based on an 
expansion of the potentials into spherical harmonics, i.e., 
into modes r ; y, m (6», </>) and r -^ l+1 ^ ^ m (6», 0), and on the 
determination of the expansion coefficients by the phys- 
ically motivated projection of the boundary conditions 
onto the modes r'Y; m (#, <j>). See also [16, Sec. 6] for 
a review of various ways to determine expansions from 
the boundary conditions in a variety of related problems. 
For nonspherical particles, our approach corresponds to 
an expansion into non-orthogonal modes and therefore is 
similar to the usage of the Raylcigh hypothesis in the the- 
ory of scattering in optics, where the scattered field at a 
perturbed interface is likewise expanded in the solutions 
of the scattered field of the unperturbed one [17]. It is 
known that such expansion methods may fail if the devia- 
tions from the ideal geometry become too large, see, e.g., 
[18] and the references therein. Nonetheless, additional 
to its simplicity and easy implementation the distinct 
advantage of our approach is its computational efficiency 
for nearly spherical particles. Thus it allows to calcu- 
late the surface modes for many realizations of randomly 
distorted nanosphcres and so to statistically characterize 
their optical responses. 

The paper is organized as follows: The numerical 
method is explained in Sec. II, and validated in Sec. Ill, 
using the cases of an ellipsoid and of a sphere with cer- 
tain shape distortions as benchmarks. In Sec. IV we give 
a statistical study of the optical response of spheres and 
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spheroids with random perturbations. 

II. THE SCHEME 

We are interested in the surface modes of sphere- 
like nanoparticles, described as some bounded domain 
li C I 3 , with boundary dil. We restrict ourselves to 
particles that are small compared to the relevant wave- 
lengths and therefore employ the quasi-static approxima- 
tion. In order to determine their surface modes we con- 
sider an excitation from infinity (described by the poten- 
tial $ext), and calculate the potential inside ($-) and 
outside ($+) the particle. These potentials fulfill the 
Laplace equation 

A<P±(x) = forx^dQ. (1) 

In addition, the boundary conditions on the surface dtl 
of the particle are 

$+(x) = ®-(x) forxedfl, (2) 
d n $+{x) = e<9„$_(£c) forzedft, (3) 

with the outward normal derivative d n and the permit- 
tivity e of the particle. The boundary condition (3) im- 
plies that the particle is surrounded by vacuum, and is 
homogeneous, isotropic, and non-magnetic; the dielectric 
properties are assumed to be local. 

There are two different methods to determine the sur- 
face modes of a nanoparticlc from (l)-(3). The first is 
to assume that the potential vanishes at infinity, i.e. to 
calculate the modes of the particle that can be present 
without an external excitation. In this case the prob- 
lem can be reformulated as an eigenvalue problem with 
a real plasmonic eigenvalue e for which a nontrivial solu- 
tion of Eqs. (l)-(3) with lim |$+(a:)| = exists [19]. 

||:c||— >oc 

In this interpretation the variable e in (3) is not regarded 
as the generally complex permittivity of the particle, but 
rather as a real eigenvalue. The second method is to 
study the system (l)-(3) with an external excitation and 
thus to regard the e in (3) as the complex permittivity 
e(uj) of the particle. The system is then solved for differ- 
ent values of the permittivity and a solution is called a 
surface mode if the field inside and around the particle 
is enhanced. If the imaginary part of the permittivity 
does not vary too much, then the enhanced fields occur 
when the real part of e is equal to a plasmonic eigenvalue. 
Thus the terms plasmonic eigenvalue and resonant value 
are closely related and will be used interchangeably. In 
general there will be a difference in the number of eigen- 
values and resonant values. While there is a infinite num- 
ber of eigenvalues, the used excitation will choose some 
of these eigenvalues and only for these an enhanced field 
will appear. 

We study the response of a particle to an external field 
and assume that the potential at infinity equals the po- 
tential of the excitation; 

lim \$ + (x)-* eBt (x)\ = 0. (4) 

||a;||— foo 



The basic idea is to expand the potential inside and out- 
side the particle into spherical modes which automati- 
cally fulfill the Laplace equation (1), 

oo / 

M*)=E E ai.mr'lHM), (5) 

;=0 m=—l 

and = + $> ex t with 

oo / 

^+o*o =E E ft,mr-< ,+i) ir(M) (6) 

1=0 m=-l 
oo I 

and $ ext (x) = E E 7;,mr<17"(M)- Here x = 

1=0 m=-l 

r(cos(0)sin(#),sin(<^)sin(0),cos(0)) and the familiar 

spherical harmonics are denoted by Y l m (9,4>) = 

P ; m (cos 9)c lm ^: For m>0, the associated Legendre poly- 

i / a \ l+m 

nomials arc P^s) = ^ — (i-^jm/a (_ j {s 2_ 1} l 

with the scale factors c] n = {-l) m sJ^^J^^,, and 

P/ n := (-l) m Pf m for m < 0. 

From (4) the coefficients 7; >m are defined by the exci- 
tation potential. Thus it remains to calculate the coeffi- 
cients ai. m and /3j )OT from (2) and (3). In order to do this 
we use the following numerical scheme. First we trun- 
cate to |Z| < N such that (N + l) 2 coefficients a^ m and 
/3/ iTO have to be calculated. To get the required 2(iV+ 1) 2 
equations we project the boundary conditions (2) and (3) 
onto the modes r l Y^ n {9 1 (f>) with degree equal to or less 
than TV. In detail, we require 

f (<S>^-^ + )r l Yr(9,<f>)dS (7) 
Jan 

= [ $Extr l Y l m (6,(t>)dS, 
Jan 

[ {ed n <S>--d n ^ + )r l Yr{6,cj))&S (8) 
Jan 

= f (d n ^> Ext )r l Y l m (e i cj ) )dS. 
Jan 

This yields a system of the form 

(Mi + eM 2 )U = M 3 G (9) 

where M U M 2 ,M 3 g C 2( - N+1 ^ x2( ' N+1 '> 2 are matrices 
which depend only on the geometry of the particle, G € 
C 2(jv+i) 2 dcpcnds only on the 7;mj and u g c 2 ( w+1 ) 2 

contains the unknown coefficients cti m and j3i m . 

For a sphere S ro of radius ro the spherical harmon- 
ics Y™ are an orthogonal (orthonormal if ro = 1) ba- 
sis of L 2 (dS ro ). Thus, (9) decouples in the case of a 
sphere (becomes block diagonal, see A for the precise 
structure) and yields (A^ + l) 2 exact solutions of (l)-(3). 
For a perturbed sphere the Y™ are no longer orthogonal, 
and the physically motivated idea of projecting onto the 
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modes r l Y™ (instead of, e.g., projecting onto the spher- 
ical harmonics Y™, which at first might appear more 
natural) is as follows: we may expect the fields to be 
localized near parts of dfl with high curvature, and for 
perturbations of spheres (of radius tq) these occur most 
naturally for parts bulging out, i.e., for r > ro. Thus, 
to minimize the error, it appears reasonable to weight 
the spherical harmonics as test functions in (7), (8) with 
r . This also complies with the folklore rule to use the 
same functions as test functions and as ansatz function. 
On the other hand, this rule is rather ambiguous here 
since we have (N + l) 2 more ansatz functions, namely 
r -{i+i)Y [ ,n , However, these tend to localize near "flat" 
parts of the perturbed sphere and are therefore less useful 
as test functions. We evaluated all three sets of test func- 
tions (r'y™k m , {Y t m )i, m , and (r^^Y™)^, against 
available exact solutions for spheroids (see Sec. Ill) and 
found that (r l Y" L ) works best while (YJ m ) and even more 
so (r - '' +1 ' Y" 1 ) yield slower convergence. 

As already pointed out in the Introduction, expansions 
like those given by Eqs. (5) and (6) are conceptually re- 
lated to the Rayleigh hypothesis, and may fail to converge 
if the deviation of the geometry considered from the ideal 
geometry is too large. Therefore it is of great importance 
to test the method against known exact solutions, and to 
control the error. As shown below, for the present prob- 
lem it turns out that already moderate N yield quite ac- 
curate results if the deviations from a sphere are not too 
large. The achievable accuracy, however, also depends on 
the quantities one wants to compute. We find that typi- 
cally N = 7, which yields Mj £ C 128x128 , is sufficient to 
calculate the resonant value of e with high accuracy. 

The generation of the matrices Mj is the most expen- 
sive part of the scheme since each entry requires the eval- 
uation of surface integrals similar to the ones in Eqs. (7) 
and (8). However, once the matrices Mj are generated, 
for any given $ est we only need to first calculate the 
coefficients 7; m and then solve some rather small linear 
system with given e. In particular, the scheme allows for 
a fast parameter scan when solving the system (9) for dif- 
ferent G, i.e., when rotating the incident field. For fixed e 
we may also define a T-matrix T = (Mi + tMq) -1 M$ to 
obtain U = TG. The calculation of the ji m is quite sim- 
ple; for instance, for a constant field in (x, y, z)*-direction 
the coefficients are -fi m = for I ^ 1 and 71,-1 = —x—iy, 
71,0 = z and 71,1 = x - iy. 

We use the GNU Scientific Library [20] for the spher- 
ical harmonics and the Cuba library [21] for calculating 
the projection integrals (7), (8). The linear system (9) is 
then solved with a standard method from LAPack [22] . 



III. COMPARISON WITH EXACT SOLUTIONS 
AND A PERTURBATION ANSATZ 

As test cases for our scheme we consider the surface 
modes of an ellipsoidal particle, for which an exact solu- 
tion exists [2] , and the case of a sphere with certain Gaus- 



sian perturbations. For the latter we compare our results 
with the results of a recently developed perturbation- 
theoretical ansatz [23]. 

A. Surface modes of an ellipsoid 

We start with a spheroid, i.e. an ellipsoid with two 
identical semi-axes. The spheroid is oriented such that 
the two identical axes are along the x- and y- axis of the 
coordinate system. Furthermore, we choose the semi- 
diameter in x- and y- direction as 1 [25] . Thus the geom- 
etry of the test case is described by one parameter R, the 
semi-axis in z-direction. As the dipole modes of a sphere 
are excitable by a constant field, and we are interested 
primarily in dipole-like modes, we use a constant incident 
field in the test cases. 

Considering a constant incident field in z-direction, the 
exact solution for the resonant value of the permittivity 
e is [2] 

e(R) = 1 - ^ (10) 

with the depolarization coefficient L(R) given by 

00 

L(R) = f / ds 1 =. (11) 

2 7 ( s + i?)V2(s + l)(s + i? 2 ) 

Numerically we determine the resonances as follows: 
After calculating the matrices Mj we solve Eq. (9) for 
different values of e[26], and define the resonant value 
as that value of e which produces the biggest dipole-like 
near field, i.e., we search the value of e which maximizes 
x /|/3i,_ 1 | 2 + |/3i, | 2 + |/3i,i| 2 . 

In Fig. 1 we compare the resonances thus obtained nu- 
merically with the exact ones, for different N and for 
two different incident fields. The results are remarkably 
good, even for small N. Indeed, for an incident field in 
z-direction and R = 1.5, for example, the exact result 
figures ase« —3.29, while the numerical procedure with 
N = 1 gives e = —3.18, and differs from the exact result 
by less than 10 -2 with N — 7. Thus, for the calculation 
of the resonant values our method gives accurate results 
even with only few spherical harmonics, and even when 
the shape deviates significantly from that of a sphere. 

On the other hand, the method is in general not well 
suited for the calculation of the fields with high accuracy 
if the deviations from the sphere become large. For illus- 
tration, we first show in Fig. 2 (a) the potentials inside 
and outside a spheroid with R = 1.4, for an incident field 
in z-direction. The linescan in (b) shows that along the 
particular line indicated in (a) the boundary conditions 
(2), (3) are reasonably well fulfilled for N = 7, while the 
jumps for N = 1 indicate that in this case N = 1 is not 
sufficient, as expected. 

From a systematic viewpoint, a more relevant basic 
check of the accuracy of the potentials is the total relative 
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FIG. 1: Comparison of the exact resonant values for a 
spheroid with the ones calculated within our approach for 
different N, and an incident field oriented along the :r-axis 
(a), and the z-axis (b). In the later Fig. 4, where we show the 
polarizability of a spheroids, the dependence of the resonant 
values on N can be seen in more detail. 

error in Eqs. (2), (3). In Fig. 3 we plot the L 2 -norms of 
these errors, again as functions of R and N, henceforth 
denoted by 

2||$+-$-|| _ 2||a n $ + -ed n $_|| 

61 ' ||* + || + ||*_||' 62 ■ llcVM + RM' 

where ||/|| = (J gn \f\ 2 AS) 1 / 2 . For R close to 1, say 0.8 < 
R < 1.3, we find that e± and are small and decrease 
monotonously in TV. However, ei,a become large rather 
quickly when R falls outside this range. Moreover, they 
then no longer decay in TV, which we attribute to a failure 
analogous to that of the Rayleigh hypothesis for such 
strongly deformed spheres. 

Similar effects can be observed for various other test 
particles: The approximation of resonant values of e is 
typically much better than e\^- Thus the performance 
of the method strongly depends on what one wants to 
compute. As a rule of thumb we find that for < 
0.2 we obtain very good approximations of the resonant 
values, and also of the polarizabilities and of the dipole 
shown below; these three observables are the 
quantities we are mainly interested in. Therefore we have 
made sure that in all calculations presented below we 
have ei t 2 < 0.1. For the solution shown in Fig. 2 with 




-2 



-2 2 

z 

FIG. 2: (Color online) (a) Calculated potentials inside and 
around a spheroid for R = 1.4 and N — 7. The imagi- 
nary part of the potential is plotted for the resonant value 
of the permittivity and the potential is normalized such that 
the maximum of the imaginary part is 1. The shape of the 
spheroid is indicated by the black ellipse. The linescans in (b) 
are taken for x — 0; the one with N = 1 is shifted downwards 
by 1 for clarity. 



N = 7 we have e2 slightly larger than 0.1, so that this 
solution would be discarded. 

Another quantity of immediate physical interest is the 
polarizability a which relates the incident field (E ext ) to 
the excited dipole moment (p e x)- In general this quan- 
tity is a tensor, but for the case of a spheroid with an 
incident field oriented along one of the principal axes the 
polarizability is a scalar [27], i.e., p ex = eoaE ext . In 
Fig. 4(a) the exact and the numerically calculated polar- 
izability of a spheroid with R = 0.8 is depicted, more 
precisely its absolute value for an incident field in z- 
direction, i.e. \a z \ cx a/I/^-iI 2 + |/3i, | 2 + |/3i,i| 2 - The 
width and the magnitude of the polarizability are ex- 
cellently reproduced by our scheme and the agreement 
between the exact solution and the numerical one is very 
good already for N > 3, and even for N = 1 the approx- 
imation already reproduces the shape of the resonance 
quite well. Fig. 4(b), which depicts the coefficients /3;. m 
for the expansion (6), shows that the dominant part of 
the field stems from / = 1 , but higher spherical harmonics 
also give non-negligible contributions. 

We have also checked our scheme against the exact so- 
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FIG. 3: (a) Z/2-norm of the boundary condition (2), and (b) 
of the condition (3), for the indicated TV. 

lution for an ellipsoid with throe different semi-diameters, 
and an incident field that is not oriented along one of the 
principal axes. Again the agreement between the exact 
solution [2] and our results is very good as long as the 
semi-diameters are not too different. 



B. Surface modes of a sphere with a perturbation 

In Sec. IV we perform a statistical analysis for the sur- 
face modes of nanoparticles with a random shape, de- 
scribed by 

r( M ) = l + S |>e X p^ 

(12) 

where s is a scaling factor, dist(#i, 4>i\ 9, fa) is the Eu- 
clidean distance between two points on the unit sphere, 
one specified by 9i and fa and the other by 9 and (f>. In 
these later studies, 9i, fa, hi and Wi will be randomly dis- 
tributed. Thus, the nanoparticles then are spheres with 
n Gaussian perturbations with height hi and width Wi. 
In order to assess our method for such cases, we first 
consider a particle with three perturbations, study how 
the resonance is shifted by varying the scaling parame- 
ter s, and compare the results with those of a perturba- 
tion ansatz for the surface modes of a nanoparticle [23]. 
For the perturbation ansatz a geometry close to the 
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FIG. 4: (Color online) Comparison of the exact and numer- 
ically calculated polarizability of a spheroid with R = 0.8 
(a). The polarizability is normalized to the one of a sphere 
with radius 1 and a dielectric function of — 2 + O.Oli. In the 
right panel (b) the coefficients of the spherical harmonics are 
presented for N = 7 and a value of e close to the resonant 
value. 

given one is needed, for which an exact solution for the 
surface modes exists. In the following that geometry is 
called the ideal geometry, and the geometry for which 
the surface modes are to be calculated is called the per- 
turbed geometry. It is assumed that the perturbation 
strength is described by a scalar parameter. In our case 
the ideal geometry is the unit sphere and the parame- 
ter that characterizes the deviation from the sphere is 
the scaling amplitude s. When making the perturbation 
ansatz the resonant values e(s) are expanded in a series 
with respect to s: 

e(s) = e(0) + se(0) + i-e(0) + • • • . (13) 

In Ref. [23] an explicit formula for e is presented against 
which we can compare our numerical results. Since the 
dipole mode of a sphere is threefold degenerate, there are 
in general three different values of e, characterizing the 
three dipole-like modes. 

The test particle with three Gaussian perturbations is 
depicted in Fig. 5(a), whereas (b) compares the results 
provided by (9) with the perturbation method. We only 
present the results for N = 7 but remark that the results 
are stable for N > 5. There is a very good agreement 
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FIG. 5: (Color online) Sketch of the used perturbed sphere 
(a), and the corresponding resonant values of e (b). The plot 
(a) refers to a scaling factor of s = 0.5; the color scale de- 
picts the distance from a point on the surface to the center 
of the underlying unit sphere. The resonant values in (b) are 
calculated with N = 7; the inset shows a blow-up near s = 0. 



between the projection and the perturbation method for 
small s. As expected there are deviations for larger scal- 
ing factors, because these are beyond the scope of the 
(first order) perturbation method. 

For s > 0.4 more than three resonances occur, notably 
one with a real part of the permittivity of about —1.6. 
By plotting the potential (or alternatively inspecting the 
expansion coefficients /3z. m ) we find that these modes are 
shifted and perturbed quadrupolc modes of the unper- 
turbed sphere, which due to the perturbations can be 
excited by a constant field. For the unperturbed sphere 
these occur at Re(e) = —1.5. Of course these modes 
are in principle present in all cases, but it depends on 
the respective geometry if these modes could be excited 
effectively by a constant field. 

To conclude the test of our scheme: Already with quite 
few spherical harmonics (typically TV < 7), and even if 
the errors ei^ in the boundary conditions are significant, 
the numerical data for the resonant values and the po- 
lariziabilitics provided by our method are stable and in 
remarkable agreement with exact solutions (if available), 
or with the results of the perturbation ansatz. 



IV. PARTICLES WITH RANDOMLY 
DISTRIBUTED GAUSSIAN PERTURBATIONS 

We now assume that we are given a set of nom- 
inally identical nanospheres which suffer from uncon- 
trolled shape imperfections induced during the fabrica- 
tion process. The task then is to characterize the optical 
properties of this set in a statistical sense. Ideally, one 
would have at least a good idea how the imperfections 
are distributed, based on an inspection of a representa- 
tive number of specimen, then generate a corresponding 
ensemble, and compute the resonances of each individ- 
ual member. Since we are not considering any particular 
case, here we simply assume that the random shape fluc- 
tuations correspond to Gaussian distortions as described 
by Eq. (12). We start with n = 4, set s = 1, and choose 
randomly distributed Bi, 0j, hi, and Wi with i = 1, . . . , 4. 
The positions of the distortions, described by the angles 
9i and are uniformly distributed on the sphere, hi is 
normally distributed with a mean of 0.2 and a standard 
deviation of 0.1; the normally distributed widths Wi have 
a mean value of fi w = 0.7 and a standard deviation of 
a w = 0.3. Clearly, perturbations with a negative width 
Wi are discarded, and we admit only perturbations with 
^f- < 2 in order to avoid sharp peaks. Therefore each 
realization is a sphere with four or less perturbations. 

In order to get significant results we generate 1000 re- 
alizations of the perturbed sphere, calculate for each re- 
alization the matrices Mi, and solve the system of equa- 
tions (9) for 100 different incident fields, i.e. we consider 
100000 cases. We choose N = 7, and require e\p < 0.1 as 
explained in Sec. Ill A, which has resulted in an ensemble 
with 42829 members. 

In Fig. 6(a) we present a histogram of the resonant 
values of the real part of the permittivity. Here we count 
a peak in the absolute value \on n \ of the polarizability 
of a particle as a resonance if \cti n \ is at least 0.1 times 
the value \ain,sphere\ of a perfect sphere. H e denotes the 
number of the resonances in the corresponding interval of 
Re(e) divided by the number of particles. In the follow- 
ing pictures H a and Hg are defined in an analogous way. 
As seen in Fig. 6(a) the shape fluctuations result in a rel- 
atively broad distribution of the resonant values around 
the unperturbed value Re(e) = —2, and the maximum is 
shifted to a slightly bigger value. 

Not only the position of the resonance is of interest, 
but also the magnitude of the induced dipole moment. 
Therefore we present in Fig. 6(b) a histogram of this 
magnitude, considering for each ensemble member only 
the resonance with the largest dipole moment. Due to the 
fact that the polarizability of a sphere scales with its vol- 
ume and a typical realization of our perturbed spheres is 
somewhat bigger than a sphere with radius 1, we normal- 
ize the polarizability to that of a sphere with the same 
volume as the perturbed one. In the large majority of 
cases the polarizability of the perturbed sphere is smaller 
than that of the unperturbed one. This is due to the fact 
that for the sphere there are no principal axes. Therefore 
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FIG. 6: Histogram of the resonant values Re(e) for the ensem- 
ble of perturbed spheres (a), and histogram of the magnitude 
of the maximum dipole moment a^™ (b). 



the polarizability is independent of the direction of the 
incident field, and for any field the induced dipole mo- 
ment points into the direction of the incident field. But 
if the symmetry of the sphere is destroyed by the per- 
turbations, there are three distinguished principal axes 
with different resonant values of the permittivity So for 
a fixed permittivity in general only one resonant condi- 
tion is matched and therefore essentially only the part 
of the incident field that points into the direction of the 
corresponding principal axis contributes significantly to 
the induced dipole moment. This results in an induced 
dipole moment which typically is smaller than that of a 
perfect sphere, and is not parallel to the incident field. 

In a plasmonic waveguide, consisting of a chain of 
metallic nanoparticles, a random angle between the in- 
duced dipole moment p and the incident field E ext influ- 
ences an efficient transport, because the near field of one 
particle should excite the neighboring particle. Therefore 
designing a plasmonic waveguide requires the knowledge 
of the spatial structure of the near field. Because a sphere 
has no principal axes, any arbitrarily small perturbation 
will pick three axes and therefore affects the orientation 
of the induced dipole. Thus, with perturbed spheres we 
may expect that there are only very few cases in which 
p is nearly parallel to E ext . To illustrate this quantita- 
tively, we present histograms of the angle between p and 
E ext in Fig. 7. In Fig. 7(a) we consider all resonances, 



whereas in Fig. 7(b) only the biggest resonance for every 
member is used. The main result is that realizations with 
p nearly parallel to E ext are indeed negligible, and the 
orientation of the dipole relatively to the external field 
is nearly random. This can be seen from the solid line 
in Fig. 7(a), proportional to sin(#), which corresponds 
to a completely random choice of the orientation of the 
dipolcs. This figure also illustrates that the only cases 
which arc suppressed when the spheres are perturbed are 
those in which p and E ext are almost perpendicular. 

As expected, when only the largest resonances are con- 
sidered the angle is typically smaller because, as already 
pointed out, in a resonant situation only that part of the 
incident field contributes that is parallel to the accord- 
ing dipole axis. Moreover, the distribution of this angle, 
shown in Fig. 7(b), is quite broad, which demonstrates 
the substantial variability in the spatial structure of the 
field around our perturbed spheres. 
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FIG. 7: Histogram of the orientation (relative to E ext ) of 
the induced dipole moment for all resonances (a), and only 
the biggest ones (b). The solid line in (a) is obtained if the 
orientation of the dipole axis is entirely random. 

Therefore, we now consider the surface modes of per- 
turbed spheroids. As spheroids already have at least one 
distinguished principal axis, we expect that perturba- 
tions of a spheroid do not have such a strong effect on 
the orientation of the induced dipoles. We employ the 
same kind of Gaussian perturbations, but now the un- 
perturbed particle is a spheroid with semi- axes R x = 1, 
Ry = 1, and R z = 1.2. Hence, the shape of a particle is 
described by 
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(sin 60 s 



cos b 
~L2 



1 -0.5 



(14) 




Q5 / dist(g. t ,0,;g,0) ' 



Again, #i and <f>i are uniformly distributed on a sphere, 
/ii are normally distributed with a mean value of 0.2 and 
a standard deviation of 0.1, and the normally distributed 
Wi have a mean value of 0.7 and a standard deviation of 
0.3. We consider 10000 realizations of the spheroid and 
calculate the response to a constant field in z-direction 
for each. Here the numerical criterion < 0.1 leaves 
us with 4323 particles. 
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FIG. 8: Location of the resonances ((a): all) and orienta- 
tion of the induced dipole moments ((b): largest only) for 
the perturbed spheroids. The vertical lines in (a), labeled 
by z and x,y, indicates the resonant values for the unper- 
turbed spheroid and an incident field along respective coordi- 
nate axes. 



the perturbed spheres. This comparison can be roughly 
quantified by the respective mean values, and the widths 
of the distributions of the angles. We focus again on the 
biggest resonances only, and determine the mean value 
9 of the angle and the interval Ie, centered around the 
mean value, which contains two-thirds of the resonances. 
For the case of the sphere this results in 6 = 1.04 and 
h = [0.65 : 1.43], whereas 6 = 0.23 and Ig = [0.09 : 0.37] 
for the case of spheroids. 



V. CONCLUSION 

We have presented a simple-to-implement and efficient 
numerical scheme for calculating important character- 
istics for surface modes of sphere-like nanoparticles in 
the quasistatic limit, based on an expansion of the "in- 
ner" and "outer" potentials into spherical harmonics. Al- 
though the spherical harmonics do not constitute an or- 
thogonal basis for particles which are not exactly spher- 
ical, and therefore encounter problems similar to those 
connected with the Rayleigh hypothesis when the de- 
viation from an exact sphere becomes too strong, they 
still remain a useful system of functions, in particular so 
when the boundary conditions are interpreted in an in- 
tegral manner, with a physically motivated choice of test 
functions. 

We have validated this scheme against exact solutions 
for ellipsoids, and against perturbation-theoretical cal- 
culations for deformed spheres. These comparisons and 
also additional intrinsic numerical tests both show that 
our method is able to yield accurate results for the res- 
onant permittivities and the polariziabilities even when 
only quite few spherical harmonics are employed, that is, 
with quite small basis sets. 

This high computational efficiency allows to perform 
statistical studies of large ensembles of randomly per- 
turbed nanoparticles. This ability is indispensable when 
designing, e.g., plasmonic waveguides from nanoparticles 
with small, but uncontrolled fabrication-induced shape 
fluctuations. On the one hand, typical effects of such im- 
perfections on the performance of these devices can be 
quantified in this manner; on the other, admissible toler- 
ances can be determined. While the specific distribution 
of shape fluctuations employed for demonstration pur- 
poses in our example given in Sec. IV may not be realis- 
tic, the key steps of such a large-scale statistical analysis 
proceed along exactly the same route as outlined there. 



In Fig. 8 we present histograms for the location of the 
resonances and for the orientation of the induced dipoles. 
Again the resonances are shifted to somewhat bigger val- 
ues of the permittivity. However, now the angle between 
the incident field and the induced dipole moment for the 
biggest resonance is much smaller, and also the distribu- 
tion of the angle is narrower than in the previous case of 
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Appendix A: Structure of the matrices 

The structure of the matrices in Eq. (9), i.e. (Mi + eM2)U = M3G, is: 



a-o, 0,0,0 


ai, 


-1,0,0 


• ajv,iv,o,o 


60,0,0,0 


61, 


-1,0,0 


6jv,jv,o,o 


cio, 0,1,-1 


0,1- 


1,1,-1 ■ 


• ajv,JV,i,-i 


60,0,1,-1 


61,- 


1,1,-1 ■ 


• 6iv,jv,i,-i 


tS0,O,Af,JV 


ai,_ 


1,N,N ■ 


• CbN,N,N,N 


bo,0,N,N 


61,- 


1,N,N ■ 


• bN,N,N,N 













do, 0,0,0 


di, 


-1,0,0 • 


■ djv,jv,o,o 










do, 0,1,-1 


di,- 


1,1,-1 ' 


• d]v,jv,i,-i 













do,o,Jv,jv 


di,. 


1,JV,JV ■ 


■ djv,jv,jv,jv 



M 2 = 



(^/,m,/' ,m' ) 



(ci, m ,i', TO ') 



where is the (AT + l) 2 x (JV + l) 2 zero matrix. The coefficients az. m ,z', m ', ■ • • , di.m,v ,m' are defined by projections 
of the boundary conditions (2), (3) onto the modes r l Y™ on the surface of the particle dfl: 



(^Lm,l' ,m' 



iLm.l'.m' 



CLmA'm' 



r(s) l+l 'Y l m (s)Y l 7 l '(s)dS , 
r(s)- {l+1 ' )+l 'Y l m (s)Y l 7 l '{s) dS . 



an 



di 



rn.V ,m' 



d n (r(s) l Y l m (s))r(sfY^'(s)dS, 

d„ (r(s)-( (+1 )y i m ( s )) /y^'(s) dS. 



on 



Thus, for a sphere with radius 1 the coefficients arc: a ; . lrM /. m / = Sa>S mm > ,h,m,v ,m' = Siv5 mm >, % m j>, m > = 

IS W S mm' 1 dl m jn,l' .m' 

= (l + l)Sii'5 mm >, and so only four block diagonals of the matrix Mi + eM 2 have nonzero en- 
tries. The vectors U and G in the system (9) are defined by the coefficients a>i im , /3/ )OT and "fi. m of the expansion of 
the potentials (5) and (6) as 



U = (aofl, ai,_i, . . . ,aN,N, Po,o, • • • , Pn,n) , G = (70,0,71,-11 ■ ■ • ,1n,n, 0, . . .,0) . 
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